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ABSTRACT 

In dense stellar clusters, binary-single and binary-binary encounters can ultimately lead to 
collisions involving two or more stars during a resonant interaction. A comprehensive survey 
of multi-star collisions would need to explore an enormous amount of parameter space, but 
here we focus on a number of representative cases involving low-mass (0.4, 0.6, and O.8M0) 
main-sequence stars. Using both Smoothed Particle Hydrodynamics (SPH) calculations and 
a much faster fluid sorting software package (mmas), we study scenarios in which a newly 
formed product from an initial collision collides with a third parent star By varying the order 
in which the parent stars collide, as well as the orbital parameters of the collision trajectories, 
we investigate how factors such as shock heating affect the chemical composition and struc- 
ture profiles of the collision product. Our simulations and models indicate that the distribution 
of most chemical elements within the final product is not significantly affected by the order 
in which the stars collide, the direction of approach of the third parent star, or the periastron 
separations of the collisions. Although the exact surface abundances of beryllium and lithium 
in the product do depend on the details of the dynamics, these elements are always severely 
depleted due to mass loss during the collisions. We find that the sizes of the products, and 
hence their collisional cross sections for subsequent encounters, can be sensitive to the or- 
der and geometry of the collisions. For the cases that we consider, the radius of the product 
formed in the first (single-single star) collision ranges anywhere from roughly 2 to 30 times 
the sum of the radii of its parent stars. The size of the final product formed in our triple-star 
collisions is more difficult to determine, but it can easily be as large or larger than a typical 
red giant. Although the vast majority of the volume in such a product contains diffuse gas that 
could be readily stripped in subsequent interactions, we nevertheless expect the collisional 
cross section of a newly formed product to be greatly enhanced over that of a thermally re- 
laxed star of the same mass. Our results also help establish that the algorithms of mmas can 
quickly reproduce the important features of our SPH models for these collisions, even when 
one of the parent stars is itself a former product. 

Key words: stars: chemically peculiar - globular clusters: general - galaxies: star clusters - 
hydrodynamics - blue stragglers - stars: interiors 



1 INTRODUCTION 
1.1 Motivations 

One exciting aspect of dense stellar systems is the simultaneous im- 
portance of three principal areas of stellar astrophysics: dynamics, 
evolution, and hydrodynamics. Many simulation codes focus on 
one of these areas and have often been lifelong works in progress. 
The first attempts at unifying these treatments into a coherent 
model to describe clusters have begun only recently. In lune 2002, 
specialists in stellar dynamics, stellar evolution, hydrodynamics, 
cluster observation, visualization, and computer science gathered at 
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the American Museum of Natural History in New York City to be- 
gin discussing a framework for Modelling DEnse STellar systems, 
without having to MODify Existing STellar codes extensively. The 
workshop-style meeting, organized by P. Hut and M. Shara, be- 
came known as MODEST- 1 (Hut et al. 2003). The second such 
meeting, MODEST-2, was organized by S. Portegies Zwart and P. 
Hut and held in December 2002 at the Anton Pannekoek Institute 
in Amsterdam (Sills et al. 2003). From MODEST-2, a set of eight 
"working groups" were established, each focusing on a different 
aspect of the MODEST endeavour.^ Attempting to integrate stellar 
dynamics, evolution, and hydrodynamics codes into one fully func- 
tional package will be challenging, largely because each area treats 
stellar properties that evolve on different time-scales. However, by 



^ See http : / / www . many body . org/modest . html. 



© 2003 RAS 



2 J. C. Lombardi, Jr., A. P. Thrall, J. S. Deneva, S. W. Fleming, and P. E. Grabowski 



combining these areas, we will be able to better model the origins, 
dynamics, evolution, and death of globular clusters, galactic nuclei, 
and other dense stellar systems. 

In this paper, our focus is on modelling hydrodynamic inter- 
actions between stars. The goal is to develop a software module 
for quickly generating collision product models, ultimately for any 
type of stellar collision, that could be incorporated into simulations 
of dense star clusters. Lombardi et al. (2002) presented an appro- 
priate formulation for treating parabolic single-single star collisions 
between low-mass main-sequence stars. Here we extend that study 
to include situations in which one of the parent stars is itself a ther- 
mally wnrelaxed collision product. Such scenarios can occur during 
binary-single or binary-binary interactions, when the time between 
collisions is much less than the thermal relaxation time-scale. 

Stellar collisions and mergers can strongly affect the overall 
energy budget of a cluster and even alter the timing of important dy- 
namical phases such as core collapse. Furthermore, hydrodynamic 
interactions are believed to produce a number of non-canonical 
objects, including blue stragglers, low-mass X-ray binaries, recy- 
cled pulsars, double neutron star systems, cataclysmic variables, 
and contact binaries. Such stars and systems are among the most 
challenging to model, but they are also among the most interest- 
ing observational markers. Blue stragglers, for example, exist on 
an extension of the main-sequence, but beyond the turnoff point. 
Blue stragglers are therefore appropriately named, as they are more 
blue than the remaining ordinary main-sequence stars, and, com- 
pared to other stars of similar mass, are straggling behind in their 
evolution. This aberration from the common path of stellar evolu- 
tion is believed to be due to mass transfer or merger in a binary 
system, or from the direct collision of two or more main-sequence 
stars (for a review, see Bailyn 1995). Predicting the nimibers, distri- 
butions, and other observable characteristics of stellar exotica will 
be essential for detailed comparisons with observations. 

1.2 Stellar dynamics and stellar evolution 

Stellar dynamics codes determine the motions of stars. The primary 
approaches to evolving clusters or galactic nuclei dynamically are 
direct N-body integrations (e.g., Portegies et al. 1997; Hurley et 
al. 2001), solving the Fokker-Planck equation (e.g., Takahashi & 
Portegies Zwart 2000), Monte Carlo approaches (e.g., Watters et 
al. 2000; Joshi, Rasio, & Portegies Zwart 2000; Joshi, Nave, & Ra- 
sio 2001; Freitag & Benz 2002; Giersz 2001; Giersz & Spurzem 
2003), and gaseous models (e.g., Spurzem & Takahashi 1995). For 
a review of the ongoing nbody effort for accurate N-body simula- 
tions, see Aarseth (1999); for a general review of cluster dynamics, 
see Meylan & Reggie (1997) and Reggie & Rut (2003). 

The most important quantities that a stellar evolution software 
module can provide to a dynamics module are the stellar masses, 
as well as the stellar radii if collisions are included. At least in prin- 
ciple, these results could come from a live (i.e., concurrent with 
the cluster dynamics) stellar evolution calculation, from fitting for- 
mulae, or from interpolation among prior calculations. Due to the 
large number of stars, it would be wasteful to expend a considerable 
amount of time for a live computation of each star's evolution. For 
the ordinary stars whose evolution is not perturbed by an event such 
as a collision, it is much more efficient and entirely appropriate to 
interpolate among, or to use analytic fitting formulae based upon, 
previously calculated evolutionary tracks. The parameter space as- 
sociated with non-canonical stars, however, is too enormous to be 
adequately covered by interpolation or fitting formulae, and it will 
ultimately be necessary to invoke a full stellar evolution calculation 



in parallel with the stellar dynamics for such stars. Although live 
stellar evolution calculations have not yet been combined with stel- 
lar dynamics codes, some parametrized codes, such as SeBa (Porte- 
gies Zwart & Verbunt 1996), SSE (Hurley, Pols, & Tout 2000) and 
BSE (Hurley, Tout, & Pols 2002), have been successfully integrated 
(Portegies et al. 1997; Portegies Zwart et al. 2001; Shara & Hurley 
2002). 

When physical collisions between stars are modelled in a clus- 
ter calculation or a scattering experiment, it is usually done using a 
method known as "sticky particles," in which a collision product is 
given a mass equal to the combined mass of the two parent stars and 
a velocity determined by momentum conservation. In a collision 
between two main-sequence stars, for example, the result would be 
modelled as a rejuvenated, thermally relaxed main-sequence star. 
This simple method is a reasonable first approximation for many 
situations. However, there are important characteristics of collision 
products that are neglected, including their rapid rotation, peculiar 
composition profiles, and enhanced sizes due to shock heating. If 
the thermal relaxation time-scale of the collision product is much 
less than the time until its next collision, then it is appropriate to 
assume the product becomes instantaneously thermally relaxed, as 
is done in the classic simulations of Quinlan & Shapiro (1990). 
This approximation becomes questionable when the collision has 
been mediated by a binary, as there is then at least one star in the 
immediate vicinity of the collision product and the likelihood of a 
subsequent collision will depend sensitively on the product's ther- 
mally wnrelaxed size. 

Binaries are subject to enhanced collision rates for two pri- 
mary reasons: (1) their collision cross section depends on the semi- 
major axis of the orbit, as opposed to the radius of a single star, and 
(2) due to mass segregation, binaries tend to be found in the core 
of the cluster, the densest and most active region. In clusters with 
a binary fraction exceeding about 20 per cent, binary-binary colli- 
sions are expected to occur more frequently than single-single and 
binary-single colUsions combined (Bacon, Sigurdsson, & Davies 

1996) . It is probably not uncommon for binary fractions to be this 
large: the inner core of NGC 6752, for example, is thought to have 
a binary fraction in the range 15-38 per cent (Rubenstein & Bailyn 

1997) . 

Binary populations can lead to complex and chaotic resonant 
interactions. These interactions tend to exchange energy between 
the binaries and the other stars in the cluster, and therefore are crit- 
ical in determining its dynamics and observable characteristics (Hut 
et al. 1992; Vesperini & Chemoff 1994). A star intruding on a bi- 
nary could, depending on parameters such as the separation of the 
binary and the velocity of the incoming star, escape to infinity, de- 
stroy the binary, form a new binary with a star from the original, or 
form a triple. The outcomes of three body encounters can be cate- 
gorized using a nomenclature based upon typical atomic processes 
(see Reggie 1975, who introduced the use of terms like "ioniza- 
tion" and "exchange," to describe resultant scenarios). See Hurley 
& Shara (2002) for an informative narration of the intimate interac- 
tions, usually involving binaries, that stars regularly undergo during 
cluster evolution. 

The STARLAB computing environment is a very useful tool for 
modelling and analyzing all types of stellar phenomena. The gen- 
eral technique for using STARLAB to determine the cross sections or 
branching ratios for the various outcomes of binary interactions is 
presented by McMillan & Hut ( 1 996) . They also highlight an exam- 
ple set of cases in which a 0.6Me star intrudes upon a binary with 
a O.8M0 primary and a O.4M0 secondary. Their assumed mass 
radius relation, M/Mq = R/Rq, is appropriate for thermally re- 
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laxed main sequence stars and is applied bothi to tlie parent stars and 
any collision products. They find that for binaries with semimajor 
axes of 0.2, O.I, 0.05, and 0.02 AU, triple-star mergers comprise 
about 1, 2, 5, and 15 per cent, respectively, of all merger events. 
As the results of the present paper will help show, we expect that 
accounting for the enhanced, thermally unrelaxed size of the first 
collision product will greatly increase these percentages as well as 
the range of semimajor axes in which triple-star mergers are signif- 
icant. 

Simulations of moderately dense galactic nuclei initially con- 
taining solar-mass main-sequence stars demonstrate that runaway 
mergers can readily produce stars with masses > lOOM©. These 
massive stars then undergo further mergers to produce seed black 
holes with masses as large as ~ 10^ M© (Quinlan & Shapiro 1990). 
This process may be responsible for massive black holes at the cen- 
tres of most galaxies, including our own. For star clusters, recent N- 
body simulations reveal that runaway mergers can lead to the cre- 
ation of central black holes within a few million years (e.g., Porte- 
gies Zwart et al. 1999; Portegies Zwart & McMillan 2002). With 
the help of Monte Carlo simulations, Rasio, Freitag, & Giirkan 
(2003) show that the runaway process will occur in a typical cluster 
with a relaxation timescale less than about 30 Myr. Observational 
evidence for a possible intermediate-mass black hole in M15 has 
been recently reported by Gerssen et al. (2002), although the data 
is more reasonably modelled with a large concentration of stellar- 
mass compact objects (Baumgardt et al. 2003; Gerssen et al. 2003). 

1.3 Stellar hydrodynamics 

Mass loss and expansion due to shock heating when two stars 
collide are examples of hydrodynamical processes that can ulti- 
mately affect the future evolution of the cluster. Mostly using the 
Smoothed Particle Hydrodynamics (SPH, see §2. 1) method, numer- 
ous scenarios of stellar collisions and mergers have been simulated 
in recent years, including colUsions between two main-sequence 
stars (Benz & Hills 1987; Lai, Rasio & Shapiro 1993; Lombardi 
et al. 1996; Ouellcttc & Pritchct 1998; Sandquist et al. 1997; Sills 
& Lombardi 1997; Sills et al. 1997, 2001, 2002; Freitag & Benz 
2003), collisions between a giant star and compact object (Rasio 
& Shapiro 1991), and common envelope systems (Rasio & Livio 
1996; Terman et al. 1994, 1995; Sandquist et al. 1998, 2000). The 
first published SPH calculations of three-body encounters were 
done by Cleary & Monaghan (1990), who performed over 100 very 
low resolution simulations and implemented a mass-radius relation 
appropriate for white dwarfs. Other three- and four-body interac- 
tion simulations include binary-binary encounters among n = 1.5 
polytropes (Goodman & Hernquist 1991) as well as neutron star 
- main- sequence binary encounters with a neutron star, main- 
sequence or white dwarf intruder (Davies, Benz & Hills 1994). See 
Rasio & Lombardi (1999) for more information concerning the use 
of SPH in stellar collisions, and see Shara (2002) for a qualitative 
overview of the progress in stellar collision research. 

If the structure and composition profiles of colliding stars were 
available (perhaps from a live stellar evolution calculation) during 
a cluster simulation, then the sticky particle method could be re- 
placed by a more detailed hydrodynamics module. SPH calcula- 
tions could then, at least in principle, be run on demand within 
this cluster simulation in order to determine the orbital trajectory 
of the product(s), as well as their structure and chemical composi- 
tion distributions. However, at least 10 ' SPH fluid particles may be 
necessary to allow an accurate treatment of the subsequent evolu- 
tion of collision products (SiUs et al. 2002). The trouble, therefore. 



is that the integration of just a single interaction could consume 
hours, days or even weeks of computing time (depending on the 
initial conditions, desired resolution, and available computational 
resources). Although the use of equal-mass particles, or the more 
accurate SPH equations of motion derived by Springel & Hernquist 
(2002), or both, could decrease the total number of particles re- 
quired, it is still currently impractical to implement a full hydrody- 
namics calculation for every close stellar encounter during a cluster 
simulation. 

One approach for incorporating strong hydrodynamic inter- 
actions and mergers into a grand simulation of a cluster, already 
successfully implemented by Freitag & Benz (2002) in the context 
of galactic nuclei, is to interpolate between the results of a set of 
previously completed SPH simulations. The SPH database of Fre- 
itag & Benz treats all types of hyperbolic collisions between main- 
sequence stars (mergers, fly-bys and cases of complete destruction), 
while also varying the parent star masses as well as the eccentric- 
ity and periastron separation of their initial orbit. The tremendous 
amount of parameter space surveyed precludes having high enough 
resolution to determine the detailed structure and composition pro- 
files of the collision products for all cases; however, critical quanti- 
ties such as mass loss and final orbital elements can be determined 
accurately. 

A second possibility is to forgo hydrodynamics simulations 
and instead model collision products by physically motivated algo- 
rithms and fitting formulae that sort the fluid from the parent stars 
(Lombardi et al. 2002). One advantage of such an approach is that 
it can handle cases in which one or both of the parent stars is itself a 
former collision product (with chemical and structural profiles that 
are substantially different than that of a standard isolated star of 
similar mass and type). 

In this paper, we use both SPH calculations and a much faster 
fluid sorting algorithm to study scenarios in which a newly formed 
collision product collides with a third parent star. By varying the 
order and orbital parameters of the collision, we investigate how 
factors such as shock heating affect the chemical composition and 
structure profiles of the collision product. Section 2 presents our 
procedures and numerical methods, both for our SPH calculations 
(§2.1) and our fluid sorting algorithm (§2.2). SPH results are pre- 
sented in §3.1, and then compared to the results of our fluid sorting 
algorithm in §3.2. In §4 we discuss our findings and possible direc- 
tions for future work. 



2 PROCEDURE 

2.1 Smoothed Particle Hydrodynamics 

One means by which we generate collision product models is with 
the parallel SPH code used in Sills et al. (2001). The original serial 

version of this code was developed by Rasio (1991), specifically 
for the study of stellar interactions such as collisions and binary 
mergers (see, e.g., Rasio & Shapiro 1991, 1992, 1994). Introduced 
by Lucy (1977) and Gingold & Monaghan (1977), SPH is a hydro- 
dynamics method that uses a smoothing kernel to calculate local 
weighted averages of thermodynamic quantities directly from La- 
grangian fluid particle positions (for a review, see Monaghan 1992). 
Each SPH particle can be thought of as a parcel of gas that traces 
the flow of the fluid, with the kernel providing each particle's spa- 
tial extent and the means by which it interacts with neighbouring 
particles. 

The SPH code solves the equations of motion of a large num- 
ber of particles moving under the influence of both hydrodynamic 
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and self-gravitational forces. All of the scenarios we investigate 
using SPH involve a O.SM© parent star, represented with 12800 
equal-mass SPH particles, and two O.GM© parent stars, each rep- 
resented with 9600 equal-mass particles. Each of our SPH particles 
therefore has a mass of 6.25 x 1O~®M0. For comparison, this par- 
ticle mass is between the masses of the central particles used in the 
AT = 3 X 10^ and Af = lO'' calculations of Sills et al. (2002), 
who used unequal-mass particles to study in detail the outer layers 
of the fluid in a collision between two O.GMq main-sequence stars. 
For our purposes, the use of equal-mass particles is more appropri- 
ate, as it allows for higher resolution in the stellar cores and does 
not waste computational resources on the ejecta. 

Local densities and hydrodynamic forces at each particle po- 
sition are calculated by appropriate summations over Nn nearest 
neighbours. The size of each particle's smoothing kernel deter- 
mines the local numerical resolution and is adjusted during each 
time step to keep Nn close to a predetermined value, 48 for the 
present calculations. Neighbour lists for each particle are recom- 
puted at every iteration using a linked-list, grid-based parallel algo- 
rithm (Hockney & Eastwood 1988). 

The hydrodynamic forces acting on each particle include an 
artificial viscosity contribution that accounts for shocks. As in Sills 
et al. (2001), we adopt the artificial viscosity form proposed by Bal- 
sara (1995), with a = f3 = 5/6, if = 0.01, and if = 10"'''. This 
form treats shocks well and has the tremendous advantage that it 
introduces only relatively small amounts of spurious shock heating 
and numerical viscosity in shear layers (Lombardi et al. 1999). 

A number of physical quantities are associated with each SPH 
particle, including its mass, position, velocity, and entropic variable 
A. Here we adopt a monatomic ideal gas equation of state, appro- 
priate for the stars in our mass range. That is, P = Ap^ , where the 
adiabatic index 7 = 5/3 with P and p being the pressure and den- 
sity, respectively. The entropic variable is closely related (but not 
equal) to specific entropy: both of these quantities are conserved in 
the absence, and increase in the presence, of shocks. 

Our code uses an FFT-based convolution method to calcu- 
late self-gravity. The fluid density is placed on a zero-padded, 3D 
grid by a cloud-in-cell method, and then convolved with a kernel 
function to obtain the gravitational potential at each point on the 
grid. Gravitational forces are calculated from the potential by fi- 
nite differencing, and then interpolated for each particle using the 
same cloud-in-cell assignment scheme. For each collision simula- 
tion in this paper, the number of grid cells is 256^. The ejecta leav- 
ing the grid interact with the enclosed mass simply as if it were a 
monopole. 

Following the same approach as in Sills et al. (2001), we begin 
by using a stellar model from the Yale Rotational Evolution Code 
(yrec) to help generate SPH models of the parent stars. We focus 
on collisions involving 0.8 and O.6M0 main-sequence stars, with 
a primordial helium abundance Y — 0.25 and metallicity Z = 
0.001. Using YREC, these stars were evolved with no rotation to 
an age of 15 Gyr, the amount of time needed for the O.8M0 star 
to reach turnoff. The total heliimi mass fractions for the 0.6 and 
O.8M0 parent stars are 0.286 and 0.395, and their radii are 0.517 
and O.955R0, respectively. See figs. 1 and 2 of Lombardi et al. 
(2002) for thermodynamic and composition profiles of the parent 
stars presented as a function of enclosed mass, as determined by 
YREC. 

To generate our SPH models, we use a Monte Carlo approach 
to distribute particles according to the desired density distribution, 
determining values of A for each SPH particle from its position. To 
minimize numerical noise, an artificial drag force is implemented. 
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Figure 1. Fractional chemical abundances (by mass) of the elements He*, 
C^^, C^^, N'^'*, O^'', Li®, Li''', and Be^ versus In A, where the entropic 
variable A = P/p^/^ is in cgs units, for our 0.6Mq (solid curve) and 
O.8M0 (dotted curve) parent stars, as determined by YREC. 



with artificial viscosity turned off, to relax each SPH parent model 
to the equilibrium configuration used to initiate the colUsion calcu- 
lations. Fourteen different chemical abundance profiles are avail- 
able from the YREC parent models to set the composition of the 
SPH particles. The abundances of an SPH particle are assigned ac- 
cording to the amount of mass enclosed by an isodensity surface 
passing through that particle in the relaxed configuration. 

Fig. 1 plots fractional chemical abundances (by mass) versus 
In A in each parent star in its relaxed SPH configuration. Note that 
the dense core of the turnoff star is at a smaller A, and its diffuse 
outer layers are at a larger A, than all of the fluid in the O.6M0 star, 
which has direct consequences for the hydrodynamics of collisions 
involving these stars. Also note that lithium and beryllium exist 
only in the outermost layers of the parent stars. 

We focus on triple-star collisions, modelling each collision 
separately and in succession. We do not consider fly-bys or graz- 
ing collisions in our SPH calculations: all of our collisions lead 
to mergers. We neglect any direct or indirect effects, including 
tidal forces, that the third star may have on the dynamics of the 
first collision. We assume that the second collision occurs before 
the first collision product thermally relaxes: a reasonable approx- 
imation since contraction to the main-sequence occurs on a ther- 
mal time-scale, lasting at least ~ lO'' yr for non-rotating products 
and as long as ~ 10^ yr for rapidly rotating products (Sills et al. 
1997, 2001), much longer than the typical time between colUsions 
in some binary-single or binary-binary interactions (but see §4.2). 

The orbital trajectory in all our collisions is taken to be 
parabolic. This is clearly not appropriate for galactic nuclei, where 
collisions are typically hyperbolic. However, in globular clusters, 
the velocity dispersion is only ~ 10 km s"'^, much less than the 600 
km s~^ escape speed from the surface of our O.8M0 turnoff star, 
and hence all single-single star collisions are essentially parabolic. 
For collisions involving binaries (even including some hard bina- 
ries), the escape speed can still be large compared to the effec- 
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Figure 2. Before initiating the second collision, we rotate the final position 
and velocity vectors of each SPH particle from the first coUision: first by an 
angle 6 around the a;-axis and then by an angle 4> around the z-axis. The 
above figure shows how a vector Q initially pointing along the z-axis is 
transformed into the new vector Q' by these rotations. 



live relative velocity at infinity. For example, con.sider two O.SM© 
tiirnoff stars in a circular orbit of radius 0.05 AU in a globular clus- 
ter. This is a hard binary, as each star moves with a velocity of 
about 60 km s~^ with respect to the center of mass of the binary, a 
speed significantly larger than the cluster velocity dispersion. Yet, 
the effective relative velocity at infinity for a collision between one 
of the binary components and an intruder would typically not be 
much more than the orbital speed, and therefore still significantly 
less than the escape speed from our turnoff star. We therefore ex- 
pect collisions between a slow intruder and a binary component to 
be close to parabolic not only for all soft binaries, but also for some 
(moderately) hard binaries 

For the first single-single star collision, the stars are initially 
non-rotating and separated by 5 Rto, where Rto = O.955R0 
is the radius of our turnoff star. The initial velocities are cal- 
culated by approximating the stars as point masses on an orbit 
with zero orbital energy and a periastron separation rp. A Carte- 
sian coordinate system is chosen such that these hypothetical point 
masses of mass Mi and M2 would reach periastron at positions 
Xi = {-iy{l-Mi/{Mi+M2))rp,yi = z^ = 0, where i = 1,2 
and i = \ refers to the more massive star. The orbital plane is cho- 
sen to be = 0. With these choices, the centre of mass resides at 
the origin. For the first collisions, the gravity grid maintains a fixed 
spatial extent from -4Rto to +4Rto along each dimension. 

For the second collision, we want to control the relative ori- 
entation of the first collision product's rotation axis and the orbital 
plane (or, equivalently, the direction of approach of the third parent 
star). To do so, we begin with the final state of the first collision and 
make two rotations to its particle positions and velocities, through 
the angles 9 and 4>. More specifically, the first rotation is clockwise 
through an angle about the x axis, while the second rotation is 
clockwise through an angle (f> about the z axis (see Fig. 2). Finally, 
the particle positions and velocities are uniformly shifted parallel 
to the x-y plane, and the third star is introduced such that the sys- 
tem's centre of mass will remain at the origin and the periastron 
positions (in the two body point mass approximation) will occur on 
the X axis. In order to allow the bulk of the fluid to remain within 
the gravity grid, the grid is extended up to a full width of IORto 
in the x and y directions for some of the second collisions. 

We use the same iterative procedure as Lombardi et al. (1996) 
to determine the bound and unbound mass. SPH structure and com- 
position profiles presented in this paper result from averaging in 
100 equally sized bins in the bound mass. Unfortunately, it is ex- 
tremely difficult to use SPH simulations to specify the equilibrium 



structure of the outermost few per cent of mass in any collision 
product. Some SPH particles, although gravitationally bound, are 
ejected so far from the system's centre of mass that it would take 
many dynamical time-scales for them to rain back onto the central 
product and settle into equilibrium. Our requirement for stopping 
an SPH calculation is that the entropic variable A, when averaged 
over isodensity surfaces, increases outward over at least the inner 
95 per cent of the bound mass in the first collision product, and at 
least 92 per cent in the second collision product. Many calculations 
are run longer in order to confirm that no rapid changes are still 
occurring in the structure and chemical composition distributions. 



2.2 Make Me A Star 

The results of parabolic collisions between low-mass main- 
sequence stars can be well explained by simple physical arguments. 
To a good approximation, the fluid from the parent stars sorts itself 
such that fluid with the lowest values of A sinks to the core of 
the collision product while the larger A fluid forms its outer lay- 
ers. Therefore, the interior structure and the chemical composition 
profiles of the collision product can be predicted accurately using 
simple algorithms, instead of hydrodynamic simulations. Based on 
these ideas, Lombardi et al. (2002) have recently created a pub- 
licly available software package, dubbed Make Me A Star (mmas).^ 
This package produces collision product models close to those of 
an SPH code in considerably less time, while still accounting for 
shock heating, mass loss, and fluid mixing. 

Sorting the shocked fluid according to its entropic variable A 
gives the A profile of the collision product as a function of the mass 
m enclosed inside an isodensity surface. In the case of the non- 
rotating products formed in head-on {Vp = 0) collisions, knowl- 
edge of the A{m) profile is sufficient to determine the pressure 
P(rn), density p(m), and radius r{m.) profiles. Using the A profile 
determined by sorting, mmas numerically integrates the equation 
of hydrostatic equilibrium with dm = Airr^pdr to determine the 
p and P profiles, which are related through p = (P/A)^^^. The 
outer boundary condition is that P = when m = Mm mas, 
where Mmmas is the desired (gravitationally bound) mass of the 
collision product. The virial theorem provides a check of the re- 
sulting profiles. This approach allows for the quick generation of 
collision product models, without hydrodynamic simulations, and 
has already been tested with single-single star collisions. 

Presented in this paper are the results from mmas for triple-star 
collisions (see §3.2). Our procedure is simple. We call the mmas 
routine twice, using the output model from the first collision as 
one of the input parent models in the second. These mmas calcu- 
lations therefore account for the differences in shock heating that 
arise from changing the order, or the periastron separations, or both, 
of the collisions. In addition to investigating all of the scenarios 
considered with the SPH code, we also use mmas to examine more 
completely how the sizes of products vary with the periastron sep- 
arations of the collisions. Furthermore, we include a O.4M0 parent 
star of radius O.357R0 , whose structure is determined by yrec un- 
der the same conditions described in §2.1. 

For an off-axis collision, knowledge of the specific angular 
momentum distribution in the collision product is necessary to de- 
termine its structure fully, which by itself is a challenging problem 
(Ostriker & Mark 1968; Clement 1978, 1979; Eriguchi & Mueller 



^ See http : / /f acuity .vassar . edu/ lombardi /mmas/. 



© 2003 RAS, MNRAS 000, 1-20 



6 J. C. Lombardi, Jr., A. P. Thrall, J. S. Deneva, S. W. Fleming, and P. E. Grabowski 



1991; Uryu & Eriguchi 1995). Although mmas outputs an approxi- 
mate specific angular momentum profile of the first collision prod- 
uct, we use only its entropic variable A(ni) profile to help initi- 
ate the second collision. That is, for one of the parent stars in the 
second collision, we always give as input to mmas the structure 
of a non-rotating star with the desired A(m) profile, a simplifica- 
tion that both eases and quickens computations. The validity of this 
approximation is supported by the SPH calculations presented in 
§3.1.3. 

Lombardi et al. (2002) implemented version 1.2 of mmas, 
while the results of this paper use version 1.6. Besides cosmetic 
changes, the primary enhancement is that the structure of the colli- 
sion product is integrated with a Fehlberg fourth-fifth order Runge- 
Kutta method. In addition, we have fine-tuned the fitted parameter 
C3 from its previous value of -1.1 to the new value -1.0, which has 
the effect of distributing shocks slightly more uniformly throughout 
the fluid. 



3 RESULTS 

Table 1 summarizes five single-single star simulations of parabolic 
collisions. The table lists: the case name; the masses Mi and M2 of 
the parent stars i = 1 and 2, respectively; the periastron separation 
rp of the initial orbit; the stellar time T/ when the calculation was 
terminated; the average radius of the isodensity surface enclosing 
90 per cent of the bound mass, as determined by SPH, Rqm.sph, 
and by mmas, R'o.9,mmas ' ^nd the final mass of the collision prod- 
uct as calculated both by SPH, Msph, and by mmas, Mm mas- 
Previous, higher resolution SPH simulations of cases e and k (see 
Sills et al. 2001; Lombardi et al. 2002) yield no significant differ- 
ences from the present calculations. For the cases in Table 1, the 
mass loss percentage ranges from about 2 to 7 per cent. Comparing 
the last two columns of this table, we see that the mass loss pre- 
scription of MMAS reproduces the SPH results for the final product 
mass to within 1 per cent in all six cases. 

Table 2 summarizes scenarios in which a collision product 
from Table 1 (referred to as the first collision product) is collided 
with a third parent star. The table shows: the case number; the name 
of the single-single collision that yielded the first collision product; 
the mass M3 of the third (i = 3) parent star; the periastron sepa- 
ration rp,2 of the second collision; the rotation angles and if) (see 
Fig. 2); the time Tf when the calculation was terminated; the ratio 
of kinetic to gravitational binding energy T/|VF| in the centre-of- 
mass frame of the final SPH collision product model; the average 
radius of the isodensity surface enclosing 90 per cent of the bound 
mass, as calculated by SPH, Ro.9,sph, and by MMAS, Ro.9,mmas', 
and the mass of the product as calculated both by SPH, Msph, and 
by MMAS, Mm MAS- 

All twenty cases presented in Table 2 involve two O.6M0 
stars and a single O.8M0 star If mass loss were neglected com- 
pletely, the mass of the final collision product would therefore sim- 
ply be 2.OM0. The SPH calculated masses range from about 1.76 
to I.9IM0, with the largest mass loss occurring for cases with suc- 
cessive head-on collisions. 

The cases in Table 2 group naturally together in a variety of 
ways. Cases 1 and 2 each involve two head-on collisions. Cases 
5 and 6 differ only in the orientation of the first collision prod- 
uct's spin axis, and an identical statement can be made for cases 7 
through 10, as well as for cases 14 through 19. Cases 2, 11, and 7 
differ only in the periastron separation of the first collision, as do 
cases 12, 13 and 14. Also, many of the cases differ only in the peri- 




„ 4x10-'" 
m 2xl0-i» 



Figure 3. Fractional chemical abundances (by mass) versus enclosed mass 
fraction m/M for a "zeroth order" collision product model generated by 
sorting the fluid from one 0.8 and two O.6M0 stars according to their A 
values, without accounting for mass loss, shock heating or fluid mixing. 
Here m is the mass enclosed within a surface of constant density and M is 
the total mass of the collision product, 2.OM0 in this model. 



astron separation for the second collision: for example, cases 7, 5, 
and 14, as well as cases 10, 6, 17, and 20. 

Even without running an SPH or MMAS calculation, one can 
generate a "zeroth order" collision product model, valid for all 
twenty cases, simply by sorting the fluid of the three parent stars by 
their A values, with A increasing from the core to the surface. In 
those regions in which more than one parent star contributes, chem- 
ical abundances can be determined by an appropriate weighted av- 
erage: the fraction of fluid with entropic variable in some range 
{A,A + A A) that originated from any one parent star is just equal 
to the fluid mass in that same range from that star divided by the 
total fluid mass in that range from all three stars. 

Fig. 3 shows the composition profiles resulting from this ex- 
ercise for the merger of two O.6M0 stars and a O.8M0 star, with 
shock heating, mass loss, and fluid mixing all completely neglected. 
The innermost 6 per cent and outermost 9 per cent of this collision 
product model consist of fluid that originated entirely in the O.8M0 
parent star, and the profiles there consequently mimic those of the 
innermost and outermost regions, respectively, of that parent. The 
profiles in the intermediate region, where all three parent profiles 
contribute, can be understood by looking back at Fig. 1 and remem- 
bering that the smaller A fluid is placed deeper in the product. For 
example, the C'^^ near m/M = 0.13 originated in the O.6M0 star, 
while the higher A, C^^-rich fluid peaked near m/M = 0.5 origi- 
nated mostly in the O.8M0 star. Although some of the composition 
profiles turn out to agree reasonably with our more precise SPH and 
MMAS calculations, the structure of the star produced by this simple 
merging procedure does not. In particular, because shock heating 
is neglected, energy is not conserved during the merger, and the re- 
sulting radius for the product is a considerable underestimate. The 
primary usefulness of the model presented in Fig. 3 is that it serves 
as a reference to help us evaluate in what ways the shock heating. 
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Table 1. Summary data of six single-single star coUisions. 



Case 


Ml 


M2 


rp 




Ro.9,SPH^ 


/?' " 

^0.9, MM AS 


MsPH 


Mmmas 




[M0] 


[M0] 


[R©] 


[hr] 


[R©] 


[R©] 


[M0] 


[Mq] 


d 


0.8 


0.6 


0.000 


6.24 


0.87 


0.88 


1.304 


1.301 


e 


0.8 


0.6 


0.368 


8.09 


1.20 


1.24 


1.369 


1.362 


j 


0.6 


0.6 


0.000 


6.24 


0.91 


0.72 


1.132 


1.121 


jk 


0.6 


0.6 


0.123 


20.9 


1.11 


0.79 


1.156 


1.149 


k 


0.6 


0.6 


0.247 


27.7 


1.52 


0.83 


1.169 


1.162 



The SPH and MMAS radii should not be directly compared. The MMAS 90 per cent radii are for 
a spherical, non-rotating star with the same A(m) profile as the collision product, while the SPH 
results account for the rotation of the product. Furthermore, the SPH radii are measured at the time 
Tf that the simulation was terminated, before aU the fluid has fallen back to the product, while the 
MMAS radii account for the shock heating this fluid wiU produce. 



Table 2. Summary data of collisions between a third star and 


a product of a first collision. 








Case 


First 


Ma 


rp,2 


e 






T/\W\ 


Iio.9,SPH 


R' " 

^0.9, MMAS 


MsPH 


Mmmas 




Product 


[M0] 


[Rq] 


n 


n 


[hr] 




[R©] 


[R©] 


[Mq] 


[Mq] 


1 


d 


0.6 


().()() 








8.78 


0.001 


1.2 


1.66 


1.765 


1.747 


2 


j 


0.8 


0.00 








15.7 


0.001 


1.3 


1.36 


1.769 


1.760 


3 


e 


0.6 


0.00 








12.0 


0.011 


1.7 


2.31 


1.799 


1.778 


4 


e 


0.6 


0.0955 








12.7 


0.059 


2.0 


3.09 


1.851 


1.825 


5 


k 


0.8 


0.198 








23.1 


0.060 


2.5 


2.05 


1.868 


1.862 


6 


k 


0.8 


0.198 


45 





23.1 


0.056 


2.4 


2.05 


1.862 


1.862 


7 


k 


0.8 


0.00 








21.0 


0.005 


1.8 


1.48 


1.794 


1.789 


8 


k 


0.8 


0.00 


90 


90 


15.7 


0.005 


1.7 


1.48 


1.788 


1.789 


9 


k 


0.8 


0.00 


45 


90 


11.5 


0.006 


1.7 


1.48 


1.793 


1.789 


10 


k 


0.8 


0.00 


45 





15.7 


0.006 


1.8 


1.48 


1.794 


1.789 


11 


jk 


0.8 


0.00 








14.8 


0.004 


1.5 


1.44 


1.792 


1.781 


12 


j 


0.8 


0.505 








23.1 


0.077 


2.4 


2.25 


1.883 


1.866 


13 


jk 


0.8 


0.505 








23.1 


0.091 


2.9 


2.48 


1.893 


1.890 


14 


k 


0.8 


0.505 








23.1 


0.092 


3.3 


2.60 


1.893 


1.902 


15 


k 


0.8 


0.505 


90 


90 


23.1 


0.081 


3.2 


2.60 


1.893 


1.902 


16 


k 


0.8 


0.505 


45 


90 


23.1 


0.088 


3.3 


2.60 


1.893 


1.902 


17 


k 


0.8 


0.505 


45 





23.1 


0.090 


3.3 


2.60 


1.895 


1.902 


18 


k 


0.8 


0.505 


90 





23.1 


0.079 


3.1 


2.60 


1.894 


1.902 


19 


k 


0.8 


0.505 


180 





23.1 


0.064 


2.7 


2.60 


1.885 


1.902 


20 


k 


0.8 


0.758 


45 





39.3 


0.091 


4.2 


2.90 


1.912 


1.917 



The same caution must be exercised when comparing the SPH and MMAS radii as with the data of Table 1. 



mass loss, and fluid mixing treated by our SPH and mmas calcula- 
tions affect the composition and structure profiles of the collision 
product. 

3.1 Smoothed Particle Hydrodynamics 

3.1.1 Varying the collision order 

Cases 1 and 2 each involve two consecutive head-on collisions be- 
tween the same three parent stars. The only variation between these 
two cases is the order in which the collisions occur: in case 1 , the 
O.8M0 parent is involved in the first collision, while in case 2 it is 
involved in the second. Fig. 4 shows the resulting cherrucal com- 
position profiles of the collision products. Because the innermost 
few per cent of the final collision products consist of low-A fluid 
that originated in the centre of the O.8M0 parent (see Fig. 5), the 
composition profiles are nearly identical there. More generally, the 



resulting profile of each chemical species is at least qualitatively 
similar throughout the products. The differences in the composi- 
tion profiles of Fig. 4 are arguably most pronounced for C^^. In the 
parent stars, this element exists in appreciable amounts only in a 
relatively thin shell, and, as shown in Fig. 1, this shell is at a higher 
value of A in the O.8M0 star than in the 0.6Mq star. Exactly where 
the C^^-rich fluid is ultimately deposited depends on the details of 
the shock heating, and hence the order in which the stars collide. 
In particular, the final C^^ profile in the case 1 product has two dis- 
tinct peaks at enclosed mass fractions of m/M ~ 0.1 and 0.5 (as in 
our zeroth order model, see Fig 3), whereas the case 2 profile has a 
single extended peak centred near m/M = 0.2. Fig. 6 displays, as a 
function of enclosed mass fraction m/M within the product, how 
each parent contributes to the overall C^^ profile. The inner peak in 
the case 1 profile is due mostly to C^^ that originated in the O.6M0 
parent stars, while the outer peak is due mostly to the higher- >1, 
C"-rich fluid from the O.SM© parent. In case 2, the O.8M0 star is 
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Figure 4. Chemical abundance profiles for the case 1 (dashed curve) and 
case 2 (dotted curve) collision products, as determined by SPH calculations. 
The chemical abundance fractions are averaged on isodensity surfaces that 
enclose a mass fraction m/M in the final collision product. 



involved in only the second collision. It therefore experiences less 
shock heating than in case 1 and more of its fluid is able to pene- 
trate to the core of the final collision product. Consequently, much 
of the C^'' from the O.GM© stars is displaced out to larger enclosed 
mass fractions, while the C^"^ from the O.SM© star is shifted in- 
ward. The net result is the single extended peak that includes C^^ 
from all three parent stars. 

The detailed differences between the composition profiles of 
the other elements in Fig. 4 can be understood by considering the 
A and composition profiles of the parent stars, along with the dis- 
tribution and amount of shock heating. For example, the core of 
the O.8M0 parent star is rich in He* and N^*, but depleted of C^^ 
(see Fig. 1). The lower shock heating to the O.SM© star in case 2 
allows more of its core to sink to the centre of the collision product. 
Consequently, the He'' and N^** levels are enhanced as compared to 
case 1, while the C^^ levels are diminished, for final enclosed mass 
fractions m/M in the range from 0.05 to 0.2. 

The mass that is ejected during the collisions comes preferen- 
tially from the outer layers of the parent stars, exactly where ele- 
ments such as Li®, IjJ , and Be'' exist. The surface abundances (by 
mass) in the final case 1 collision product for these three elements 
are approximately 3 x 10"^°, 4 x 10"®, and 2 x 10"^", which 
is, respectively, about 30, 6, and 3 times less than at the surface of 
the O.8M0 parent star. In case 2, the surface layers are compara- 
bly depleted in these elements: the corresponding abundances are 
2 X 10"^°, 6 X IQ-^, and 2 x 10"^°. 

The bottom three panes in Fig. 4 show that the distribution 
of Li", Li^, and Be'* does differ somewhat between cases 1 and 2. 
Because the CSM© star suffers less shock heating in case 2 than 
in case 1, it loses less of its mass as ejecta and, consequently, can 
contribute more Li®, Li'^, and Be^ to the outermost layers of the 
final collision product. Furthermore, when the fluid containing Li^ 
and Be'* from the O.6M0 stars is involved in both collisions (case 
2), it is shocked more and ultimately either ejected or deposited in 



Figure 5. Fractional contribution fi to the mass of the collision product 
from each parent star as a function of enclosed mass fraction m/M within 
the product for case 1 (top) and case 2 (bottom), as determined by SPH cal- 
culations. Each of these cases involves head-on (r-p = 7-^,2 = 0) collisions 
among one 0.8 and two O.6M0 stars; however in case 1 the O.SMq star is 
part of the initial collision, whereas in the case 2 scenario it is part of the 
second coUision. Different line types are used for each parent star: i = 1 
(solid curve), i = 2 (dotted curve), and i = 3 (dashed curve). Parents with 
an index i = 1 or 2 are involved in the first collision, while i = '?> refers 
to the third parent star from the second collision. The contribution profile 
from the O.8M0 parent is labelled in each case. 



the outer ~10 per cent of the product. However in case 1, the one 
O.6M0 star that is involved in only a single collision can deposit 
its Li^ and Be'' of comparatively \ow-A deeper in the product, re- 
sulting in flattened profiles extending further into the interior. 

Fig. 5 shows, as a function of m/M, the fractional contribu- 
tion to the final product's mass from each of the three parent stars 
for cases 1 and 2, as determined by SPH calculations. In each case, 
the innermost few per cent of the final collision product consists of 
low-A fluid that originated in the centre of the O.SM© parent. Be- 
cause the first collision in cases 1 and 2 is head-on (rp = 0), fluid 
from the first two parent stars is not distributed axisymmetrically in 
the first collision product (the composition distribution is therefore 
not axisymmetric, even though the structure of the first product is). 
In case 2, the O.SM© star strikes the first collision product on the 
side with fluid from the first (i = 1) O.6M0 parent. Fluid from the 
first O.6M0 parent is therefore heated more than fluid from the sec- 
ond, and the former is buoyed out to larger enclosed mass fractions 
in the final product. In off-axis collisions, rotation induces shear 
mixing, so that if two identical stars are involved in the first col- 
lision, they contribute essentially equally within the final product: 

/l = /2. 

The profiles of Fig. 4 and Fig. 6 demonstrate that the order in 
which the stars collide can influence shock heating enough to af- 
fect, at least slightly, the chemical composition distribution within 
the final collision product. While the difference in resulting chemi- 
cal composition profiles is small, Fig. 7 shows that the difference in 
the structure of the collision product would be completely negligi- 
ble for most purposes. Although changing the order of these head- 
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Figure 6. Fractional abundance of C^"* versus the enclosed mass fraction 
m/M in the final collision products of case 1 (dashed curve) and case 2 
(dotted curve), as determined by SPH calculations. The top pane shows 
the total C^^ abundance. The second pane shows the contribution from the 
0.8Mq parent, and the bottom two frames show the contribution from the 
0.6Mq parents, so that the curves in the bottom three panes add up to give 
the overall profile shown in the top pane. 



Figure 7. Structural profiles for the case 1 (dashed curve) and case 2 (dotted 
curve) collision products, as determined by SPH calculations: the enclosed 
mass fraction m/M, the natural logarithm of the average entropic variable 
A, and the base 10 logarithm of the average density p, are aH plotted as a 
function of the average distance r from the centre of the coUision product 
to an isodensity siuface. Units are cgs. 



on collisions affects how the shock heating is distributed (and hence 
where any particular fluid element settles), it does not greatly affect 
the overall amount of heating that occurs nor the amount of mass 
that is ejected. At least for low Mach number collisions (as with 
parabolic collisions) that are nearly head-on (so that the merger 
occurs quickly), shock heating can be thought of as a mild pertur- 
bation; consequently, the final A profile, and hence the structure of 
the final product, is not sensitive to the collision order in such cases. 



3.1.2 Varying the direction of approach of the third parent 

We now investigate how the direction of approach of the third star 
toward the first collision product affects the final collision product. 
One might wonder, for example, whether an impact in the first col- 
lision product's equatorial plane (^ = or 180°) would yield a 
qualitatively different result than if the impact had instead occurred 
on the rotation axis. Cases 5 and 6, cases 7 to 10, and cases 14 to 19 
can all be used to examine such effects, as the cases within each set 
differ only in the angles 6 and 0, by which the first product is ro- 
tated (see Fig. 2). We find that while the spin of the final product is 
of course sensitive to such variations (e.g., see the T/IM^I colurrm 
of Table 2), the composition profiles are nearly unaffected. 

Fig. 8 shows the chemical abundance profiles of the collision 
product resulting in four cases in which the angle of approach of the 
third star is varied (cases 14, 15, 17, and 19). Each of these cases 
involves an off-axis collision between a case k collision product and 
a O.8M0 star. In case 14, the first collision product's spin vector is 
parallel to the orbital angular momentum of the second collision. In 
the other cases, the case k collision product is tilted various ways 
according to the values of 6 and 4> listed in Table 2. In case 19, for 




a. 2x10-1° 

m 10-'° 




Figure 8. Chemical abundance profiles for the final collision product of 
cases 14 (solid curve), 15 (dotted curve), 17 (long dashed curve), and 19 
(short dashed curve), as determined by SPH calculations. In these scenarios 
the same rotating first collision product collides off-axis with a O.SMq star, 
with a different orientation of the first collision product's rotation axis in 
each case. 
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Figure 9. Fractional contribution f j, from tlie tliird parent star as a func- 
tion of enclosed mass fraction m/M within the collision product of cases 
14 (solid curve), 15 (dotted curve), 17 (long dashed curve), and 19 (short 
dashed curve), as determined by SPH calculations. In these scenarios, the 
first two parent stars are both O.6M0; the third parent star is 0.8Mq and 
approaches from a different angle 6 relative to the rotation of the first colU- 
sion product in each case. The fractional contribution from each of the first 
two parent stars is essentially equal and can therefore be determined easily 
from the /a curve: /i = /2 = (1 - /3)/2. 



example, the case k collision product is flipped over 180° so that it 
rotates in an opposite direction to that of the case 14 rotation. 

In case 14, the fluid of the first product is rotating with the third 
star's motion as it impacts {9 = 0°). Consequently, the merger pro- 
cess is relatively gentle. For larger 0, the relative impact velocity is 
larger and the merger is somewhat more violent. Cases 14, 17, 15, 
and 19 have 9 values of 0, 45, 90, and 180°, respectively; as 6 in- 
creases, slightly less fluid from the O.SM© star can sink down into 
the core of the final collision product, and the C^^ profile rises at 
a slightly smaller enclosed mass fraction m/M (see Fig. 8). Nev- 
ertheless, as shown in Fig. 9, the contribution of each parent star 
to the product varies very little from case to case. Consequently, 
the chemical profiles in the collision products also vary little as 9 
is changed. Indeed, the He", C^^, N^"*, and O"' profiles in Fig. 8 
all look remarkably similar to the corresponding profiles in Fig. 3 
for our simple, zeroth order model. However, the C^^ profile has a 
single broad peak, for the same reasons as in case 2. Furthermore, 
because of mass loss, the beryllium and lithium surface abundances 
are found to be greatly less than our zeroth order model (that ne- 
glects mass loss) would indicate. 

The structure of the final collision product (see Fig. 10) can 
be affected by the direction of approach for two primary reasons. 
Firstly, having larger relative velocity at impact leads to larger 
shock heating. Notice, for example, how the case 19 product has 
the largest A values in Fig. 10. Secondly, having less angular mo- 
mentum in the system leads to a more compact product. For ex- 
ample, the case 19 product has the largest enclosed mass fraction 
for almost any average radius r, despite the additional shock heat- 
ing undergone in this case. Furthermore, by comparing the final 




Figure 10. Structural profiles for the final collision product of cases 14 
(solid curve), 15 (dotted curve), 17 (long dashed curve), and 19 (short 
dashed curve), as determined by SPH calculations. The particular quanti- 
ties plotted are as in Fig. 7. 

masses Msph listed in Table 2 for the products of cases 5 and 6, 
of cases 7 through 10, and of cases 14 through 19, one can see that 
the amount of mass ejected is only very weakly dependent upon 
the direction of approach of the third parent star, varying by about 
O.OIM0 or less within each of these sets of cases. 

3.1.3 Varying the periastron separations of the collisions 

We now investigate the effects that the periastron separation of the 
first collision has on the final collision product. Cases 12, 13, and 
14 all involve off-axis collisions with first colUsion products that 
are created from the same O.GM© parent stars but with different 
periastron separations (cases j, jk, and k, respectively), and hence 
different inherited angular momenta. Fig. 11 plots the fractional 
contribution of the third parent star within the merger product. In all 
cases, the low-A core of the O.SM© star is able to sink to the core 
of the final product. However, as the periastron separation of the 
first collision is increased, the two O.6M0 parent stars experience 
more shock heating, and the O.SMq parent is able to have more 
fluid penetrate down to the depths near m/M ^ 0.15. 

Fig. 12 presents chemical composition profiles of the final 
collision product for these three cases. These profiles demonstrate 
that the angular momentum of the first collision product has only 
a small effect on the final collision product. As expected from 
Fig. 11, the profiles of the three cases are essentially identical in 
the innermost 5 per cent of the bound mass, because only the core 
of the 0.8 Mq parent contributes there. The abundance profiles of 
each chemical species are at least qualitatively, and usually quanti- 
tatively, similar throughout the products. The variations that do ex- 
ist can be understood in terms of the different shock heating dtiring 
the first collisions. Because the amount of shock heating increases 
with periastron separation, the case j, jk, and k products have in- 
creasingly larger values of A at almost any enclosed mass fraction 
(this trend is not immediately evident in Fig. 13 only because In A 
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Figure 11. Fractional contribution of tlie third parent star, as determined 
by SPH calculations, within the final collision product for three scenarios 
that differ only in the periastron separation rp of the first collision: cases 12 
(solid curve), 13 (dotted curve), and 14 (dashed curve). In all three cases, 
the first two parent stars are O.6M0, while the third parent is O.SM©. The 
O.SMq parent penetrates into the product the least in case 12, because of 
the relatively small amount of shock heating suffered by the first product 
during the first coUision in this case. 




Figure 12. Chemical abundance profiles, as determined by SPH calcula- 
tions, for the final collision product for three scenarios that differ only in 
the periastron separation of the first coUision: cases 12 (solid curve), 13 
(dotted curve), and 14 (dashed curve). 



Figure 13. Structural profiles, as determined by SPH calculations, for the 
final collision product for three scenarios that differ only in the periastron 
separation rp of the first collision: cases 12 (solid curve), 13 (dotted curve), 
and 14 (dashed curve). The particular quantities plotted are as in Fig. 7. 



is being plotted versus radius and not enclosed mass). The fluid 
from the O.SM© star is therefore able to penetrate the case j prod- 
uct the least, the case jk product a little more, and the case k product 
even more still (see Fig. 11). Consequently the rise in C^^ and C^^ 
abundance is pushed out to increasingly larger enclosed mass frac- 
tions rn/M in Fig. 12 as one considers cases 12, 13, and 14, in that 
order. In case 12 and arguably case 13, the cases with the lesser 
amounts of shock heating, traces of two separate peaks are evident 
in the C" profile. As in our zeroth order model (see Fig. 3), the 
inner peak is due mostly to C^^ from the O.6M0 stars while the 
outer peak is mostly due to C^^ from the O.8M0 star. 

Fig. 13 shows that the structure of the bulk of the fluid in the 
final collision product is not significantly affected by the periastron 
separation of the first collision, and hence the spin of the first colli- 
sion product. There is, nevertheless, a visible trend for the enclosed 
mass fraction at a given average radius to decrease for products 
with more spin. For example, the isodensity surface with an av- 
erage radius of 3.2Rq encloses about 94 per cent of the case 12 
product, about 92 per cent for the case 13 product, and only about 
90 per cent of the case 14 product. Such a trend is expected, simply 
because of expansion due to rotational support. 

Cases 10, 17, and 20 can be used to investigate the effects that 
the periastron separation of the second collision has on the profiles 
of the final product. Cases 10, 17, and 20 involve collisions between 
a case k product and a O.SM© star, with periastron separations for 
the second collision of rp,2 = 0, 0.505, and O.TSSR©, respectively, 
corresponding to a number of passages or interactions rip = 1, 2, 
and 3, again respectively. See the discussion of fig. 6 in Lombard! 
et al. (1996) for the details of how rip is determined. 

Fig. 14 reveals the way in which the O.SM© parent contributes 
to the final product in each of these three cases. As usual, the low- 
A core of the O.SMo star sinks to the core of the collision product. 
As the periastron separation of the second collision is increased, 
the resulting colUsion products tend toward larger mass-averaged 
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Figure 14. Fractional contribution of the tliird parent star, as determined by 
SPH calculations, within the final collision product for three scenarios that 
differ only in the periastron separation rp,2 of the second collision: cases 
10 (solid curve), 17 (dotted curve), and 20 (dashed curve). In aH three cases, 
the first two parent stars are O.6M0, while the third parent is O.SMq. The 
first coUision product results from case k. 



Figure 15. Chemical abundance profiles, as determined by SPH calcula- 
tions, for the final collision product for three scenarios that differ only in the 
periastron separation r-p,2 of the second coUision: cases 10 (solid curve), 17 
(dotted curve), and 20 (dashed curve). 



values of A. Fluid from the the O.8M0 star therefore can penetrate 
the case 10 product the most, the case 17 product a little less, and 
the case 20 product even less still. 

Fig. 15 shows that the composition profiles in these three cases 
are essentially identical in the innermost few per cent of the col- 
lision products. Indeed, the abundance profiles are again at least 
qualitatively, and usually quantitatively, similar throughout the en- 
tire product. As before, slight variations do result from having dif- 
ferent distributions of shock heating. In particular, the rise in C^^ 
and C^^ abundance is drawn in to smaller enclosed mass fractions 
m/M in Fig. 15 as one examines cases 10, 17, and 20, in that or- 
der. Fig. 16 reveals the differences in the final product structure for 
these three cases. The top pane shows that the mass distribution 
of the final product is affected by the periastron separation of the 
second collision in a way that is simple to understand: increasing 
the second periastron separation increases both shock heating and 
rotation, and so a given radius encloses less mass. 

3.2 Fluid sorting witli mmas 
3.2.1 Comparison with SPH results 

In §3.1.2 we found that the direction of approach of the third star 
only weakly affects the profiles and mass of the final product. We 
therefore do not account for the angles 9 and </> of the second col- 
lision when applying our fluid sorting package mmas. As a result, 
the product model that mmas generates is identical within each of 
the following sets: cases 5 and 6, cases 7 through 10, and cases 14 
through 19. In all twenty cases presented in Table 2, the final prod- 
uct mass given by mmas agree with those from SPH to within 1.5 
per cent. 

Fig. 17 compares the chemical composition profiles of final 
collision products, as determined by both mmas and SPH models. 




Figure 16. Structural profiles, as determined by SPH calculations, for the 
final collision product for three scenarios that differ only in the periastron 
separation rp,2 of the second colUsion: cases 10 (sohd curve), 17 (dotted 
curve), and 20 (dashed curve). 



for two scenarios (cases 1 and 2) in which each collision is head-on 
(rp = rp,2 = 0). These cases involve the same three parent stars; 
however the order in which the stars collide is varied. The mmas 
abundance profiles maintain the same qualitative shape as those of 
the SPH data for almost all of the elements. One possibly important 
difference is that the mmas package slightly over-mixes the core in 
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Figure 17. Chemical abundance fractions by mass vs. tlie enclosed mass fraction m/M in the final collision product for (a) case 1 and (b) case 2. Results both 
from an SPH simulation (dotted curve) and from the MMAS software package (solid curve) are shown. In each case, both colUsions are head-on. 





Figure 18. Structure profiles as a function of radius r in the final collision product for (a) case 1 and (b) case 2. Results both from an SPH simulation (dotted 
curve) and from the MMAS software package (solid curve) are shown. In each case, the final coUision product is non-rotating. 



case 1, and, consequently, the central helium abundance is not quite 
as high as in the SPH calculation. Another noteworthy difference is 
that the Li^ profile, especially in case 1, is not well represented near 
the surface. Because Li® exists in an even thinner shell at the sur- 
face of the 0.8M© star than does Ia' and Be®, its abtmdance profile 
in the product is particularly sensitive to the mass loss distribution 
during the collisions. Note that mmas does correctly predict that 
most Li^ is ejected during the collisions. Furthermore, the abun- 
dance profiles generated by mmas much more closely resemble the 
SPH results than our zeroth order model does (see Fig. 3), indicat- 



ing that MMAS is capturing the important effects of mass loss and 
shock heating. 

In scenarios such as cases 1 and 2 for which the final product 
is non-rotating, it is straightforward to obtain the enclosed mass m 
and density p profiles from the A profile by integrating the equation 
of hydrostatic equiUbrium (see §2.2). Fig. 18 shows the resulting 
structure of the final collision products. The kink in the A profile 
a little inside r = O.IR0 marks the boundary within which fluid 
from only the O.8M0 star contributes, and MMAS reproduces this 
feature quite well. The central density of the SPH model is slightly 
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Figure 19. Chemical abundance fraction by mass versus the enclosed mass 
fraction m/M m the final collision product for both the MMAS (solid curve) 
and the SPH (dotted curve) data of the case 4 collision product. 



less than that of the mmas model, mostly due to how density is cal- 
culated as a smoothed average in SPH. Despite this difference, the 
overall structure of the collision product is extremely well repro- 
duced by MMAS. 

Fig. 19 compares the chemical composition profiles for SPH 
and MMAS data for case 4, a situation in which both collisions 
are off-axis. The most noticeable discrepancies are that mmas 
again slightly over-rrrixes the core and underestimates the surface 
Li® abundance. Nevertheless, the chemical abundance profiles pro- 
duced by the mmas package and the SPH code are extremely sim- 
ilar. For example, mmas correctly reproduces the C^^ abundance, 
with three peaks each corresponding to a different parent star. The 
inner peak is due to the low- A fluid from the O.6M0 star involved 
in only the second collision, the middle peak represents fluid from 
the other 0.6Mq star, and the outer peak represents high-A fluid 
from the 0.8Mq parent (see Fig. 20). Note that this feature is re- 
producible by mmas only because it accounts for shock heating 
in each collision (compare to our zeroth order model of Fig. 3, in 
which there are only two peaks in the C^"' profile). 

In many of the MMAS models, small kinks, or discontinuities, 
are evident in some of the abundance profiles: such features mark 
locations outside of which an additional parent star either starts or 
stops contributing. For example, in the case 5 and 6 collision prod- 
uct, a kink exists in the C^'^ and C^'^ profiles near rn/M = 0.08 
(see Fig. 21). As is evident from Fig. 22, fluid inside of the 
m/M « 0.08 shell originated solely in the O.8M0 parent star. 
In the range m/M > 0.08, all three parent stars contribute. The 
smoothing that is inherent to the SPH scheme makes it difficult to 
resolve such features with our hydrodynamics code. It is possible 
that similarly abrupt changes in abundance could occur in nature 
within real collision products. 

By comparing the SPH data within Fig. 21, as well as within 
Fig. 22, wc also see an example of the trend discussed in §3.1.2. 
Namely, the direction of the first collision product's rotation axis 
(or equivalently the direction of approach of the third parent) has 



Figure 20. Fractional abundance of C^"^ versus the enclosed mass fraction 
m/M in the final case 4 collision product, as determined both by MMAS 
(soUd curve) and by SPH (dotted curve). The top pane shows the total C'^^ 
abundance, while the bottom three panes show the contributions from each 
individual parent star: Cj^ is the contribution from the 0.8Mq parent, 
is the contribution from the O.6M0 parent in the first collision, and Cg^ is 
the contribution from the 0.6Mq parent in the second collision. 
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Figure 21. Chemical abundance fraction by mass versus the enclosed mass 
fraction m/M in the final coUision product for both the MMAS (solid curve) 
and the SPH data of the case 5 (dotted curve) and case 6 (dashed curve) 
collision products. Note that our MMAS results do not distinguish between 
cases 5 and 6, as they differ only in the orientation of the first product's spin. 
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Figure 22. Fractional contribution fi of each parent star i versus the en- 
closed mass fraction m/M in the case 5 and 6 collision products, as de- 
termined both by MMAS (solid curve) and by SPH. The dotted curve corre- 
sponds to the case 5 SPH results, while the dashed curve gives the case 6 
SPH results. The same MMAS model is valid for both cases, as they differ 
only in the direction of the first product's rotation axis. The i = 1 and 2 
parents are O.6M0, while the i = 3 parent is O.SM©. 

little effect on the final collision product. Indeed, when using mmas, 
our approach is to neglect completely the rotation of the first colli- 
sion product, which is why the same mmas model applies to both 

cases 5 and 6. 

Fig. 23 plots the entropic variable A versus the enclosed mass 
fraction for the final collision products of cases 18, 19, and 20, as 
determined both by SPH and mmas. Cases 18 and 19 differ only 
in the direction of approach of the third star, and we again see that 
this variation has little effect on the SPH results. The kink in all 
of the profiles slightly inside rn/M = 0.1 marks the boundary 
within which fluid from only the O.SMq star contributes. MMAS 
again reproduces this feature quite well, mmas does underestimate 
the shock heating to the core and hence the central value of A, 
although some of this discrepancy is due to the spurious heating 
evident in longer SPH simulations (Lombard! et al. 1999). Never- 
theless, it is likely that this difference between the mmas and SPH 
models would last only as a transient during the thermal relaxation 
in a stellar evolution calculation. It is also worth noting that, while 
the SPH calculations need to be terminated before all of the bound 
fluid can settle into equiHbrium (see §2.1), the mmas A profile does 
steadily increase outward throughout the entire product. 

3.2.2 Sizes of collision products 

Unfortunately it is a difficult task to determine the overall size of 
a colUsion product, either with SPH simulations or with a package 
like mmas. Whenever there is any mass loss in an SPH simulation, 
there will also be SPH particles that are nearly unbound and, in 
practice, still moving away from the product when the simulation 
is terminated. These particles would ultimately form the outermost 
layers of the collision product, but it would take an utterly unfea- 



sible amount of time to wait for them to come back and settle into 
equilibrium. 

The entropic variable A profile produced by mmas seems quite 
reasonable, both because it increases aU the way out to the surface 
and because the SPH results tend to approach its form as more of 
the fluid settles into equilibrium. However, there are no simulation 
data to compare against for the very outermost layers of a product 
and so the exact form of the profile there is difficult to validate. Not 
surprisingly, the radius of the collision product is rather sensitive 
to the A profile. For example, simply by changing the parameter 
C3 from -1.0 to the still very reasonable value of -1.1, which tends 
to distribute slightly more shock heating to the outer layers (see 
Lombardi et al. 2002), the radii of our mmas final collision prod- 
uct models for case 1 and for case 2 increase by about a factor of 
2. Despite such uncertainties, it is still interesting to get a crude 
estimate of the sizes of the collision products immediately from 
MMAS. In making these estimates, we do not account for the expan- 
sion due to rotation, but instead simply integrate the equation of 
hydrostatic equilibrium for a non-rotating star with the same A pro- 
file, using the outer boundary condition that the pressure vanishes. 
The radii calculated therefore represent the sizes that the products 
would have if some mechanism were to brake their rotation without 
disturbing their A profiles. 

Fig. 24 plots the radii at various enclosed mass fractions for 
products generated in single-single star collisions involving 0.4, 0.6 
and O.8M0 parent stars, as determined by mmas. These radii are 
plotted against the normalized periastron separation Vp / {Ri + R2) , 
which we allow to exceed unity slightly to account for bulges in the 
parent stars. The general trend is that as the periastron separation 
increases, the collisions are more long-lived, there is more shock 
heating, and the radii of the collision products increase. Because 
the fluid in the deep interior of the product is largely shielded from 
shocks, the A profile there, and hence the radius r profile, are not 
too strongly dependent on the periastron separation of the collision. 
As a result, the radii versus periastron separation curves of Fig. 24 
become closer to horizontal as one looks to smaller enclosed mass 
fraction. For the cases examined in Fig. 24, the full (100 per cent 
enclosed mass) radius of the collision product is always at least 
about twice the sum of the radii of the parent stars, and often even 
much larger than this. For example, if two O.8M0 stars suffer a 
grazing (rp Ri + R2) collision, the collision product then has a 
full radius of about 4OR0, about 20 times larger than the sum Ri + 
R2 of the parent star radii. We therefore expect that the collisional 
cross-section of these first products will be significantly enhanced 
over that of their thermally relaxed counterparts. 

Fig. 25 is similar to Fig. 24, but for triple-star collisions. We 
use different line types to represent various normalized periastron 
separations rp/{R\ + R2) for the first collision, and along the hor- 
izontal axis we vary the normalized periastron of the second col- 
lision. The curves give the radii at three different enclosed mass 
fractions. For each of the six r-p values in a frame of Fig. 25, we 
performed a nested loop over 45 equally spaced values of the nor- 
malized periastron separation for the second collision, from to 
1.1. Therefore, mmas treated 270 different triple-star collisions (in 
a few minutes on a Pentium IV workstation) for each of the four 
plots. Note that there is a general trend for the radius of the colli- 
sion product to increase as the first periastron separation increases, 
as expected; this effect is mild for the 50 per cent enclosed mass ra- 
dius, and rather dramatic for the full radius. Even more significant 
is the second periastron separation, with grazing second collisions 
resulting in products that are substantially larger than those from 
head-on collisions: the shock heating suffered by the already dif- 
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Figure 23. Entropic variable A as a function of enclosed mass fraction m/M in representative final collision products. In frame (a), the same MMAS model 
(solid curve) is compared against our SPH models for cases 18 (dotted curve) and 19 (dashed curve): our implementation of MMAS does not distinguish 
between these cases, as they differ only in the orientation of the first collision product's rotation axis. In frame (b), the MMAS (solid curve) and SPH (dashed 
curve) results are compared for case 20. The MMAS profiles have the same qualitative form as the SPH results, except in the outer ~ 10 per cent of the bound 
mass where the SPH models have not settled into equilibrium. 




Figure 24. As a function of the normalized periastron separation rp/(Ri + R2), each pane shows the radius R' that encloses, from the top curve to the 
bottom one, 100, 99, 95, 86. and 50 per cent of the total bound mass of the collision product, on a logarithmic scale and as determined by MMAS. The scale 
on the left gives the radius in solar units, while the scale on the right normahzes the radius to the sum of the parent star radii. The combination of parent stars 
considered are (a) 0.8 and O.8M0, (b) 0.6 and 0.6M©, (c) 0.4 and O.4M0, (d) 0.8 and 0.6Mq, (e) 0.8 and O.4M0, and (f) 0.6 and O.4M0. 
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Figure 25. We plot various radii R' of tlie final collision product (as determined by MMAS) as a function of the periastron separation of the second collision, 
normalized to sum of the radii of the colliding stars. The radius ^?'i02 ^^^^ collision product is the 100 per cent radius from Fig. 24 corresponding to 

various normalized periastron separations for the first collision: r-p i/ (Ri + R2) = (solid curve), 0.2 (dotted curve), 0.4 (short dashed curve), 0.6 (long 
dashed curve), 0.8 (dot - short dashed curve), and 1.0 (dot - long dashed curve). There are three curves of each line type: the bottom one corresponds to the 
radius enclosing 50 per cent of the final bound mass, the middle one coiTesponds to the 95 per cent radius, and the top one to the 100 per cent radius. 



fuse outer layers of the first collision product is severe when multi- 
ple pericentre passages occur before merger. Once rp,2 grows large 
enough for the third star's initial impact to be outside of the first 
product's core, so that more than one pericentre passage would re- 
sult before merger, then the shock heating is no longer as sensi- 
tive to rp,2 and the full radius surfaces in Fig. 25 tend to plateau. 
How strongly the full radius varies with the rp,2 therefore depends 
on the mass distribution within the first product. For first products 
with a more uniform density, such as in the product of two O.4M0 
stars, the final product size increases more gradually and consis- 
tently with rp,2- 



As the mass of any one of the three parent stars is increased, 
the trend is for the radius of the collision product to increase as 
well. For example, Fig. 25(a) shows that for collisions in which two 
0.6Mq stars collide and then a O.SM© collides with the first prod- 
uct, the final collision product radius does not exceed a few times 
IO^Rq. If one of the O.6M0 stars is substituted with a O.8M0 
star, then the final radius can be as large as about lO^'R© [see 
Fig. 25(b)]. This extreme size is due to the phenomenally diffuse 
outer layers of the product: the average density of such a star is 
only ~ 10^^* g cm^'^. The noise visible on some of the full ra- 
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dius curves is due to approaching ttie limiting numerical precision 
during the structure integration in these diffuse regions. 

From Fig. 25 we see that the radius that encloses 95 per cent 
of the total mass, while still large, is often orders of magnitude 
smaller than the full radius of the final product. Because of the low 
densities involved, the full radius calculated is rather sensitive to 
the details of the shock heating during the collision. Changing the 
MMAS parameter cz from -1.0 to -1.1, for example, can increase 
the full radius by a factor of a few, although the radius enclosing 
50 per cent of the total mass does not change by more than a few 
per cent. Nevertheless, any reasonable form and amount of shock 
heating yields products that are significantly larger than a thermally 
relaxed star with the same mass and composition. 

Colliding the same three parent stars in a different order does 
not drastically affect the mass of the final product, although it does 
significantly affect its size. Consider, for example, frames (c) and 
(d) of Fig. 25. If two 0.4Mo stars collide and then the resulting 
product collides with a O.SMq star, the final product t5^ically has 
a radius of order 10 — lOOR©, but if the O.SMo star is switched 
into the first collision instead, the final radius is usually in the range 
from ~ 100 — 10^ R© . The primary reason for this difference is that 
a collision between the 0.4 and O.8M0 stars yields a product with 
especially diffuse outer layers, and, as a result, is subject to a larger 
number of passages and hence more shock heating during a second 
collision. 



4 DISCUSSION 

4.1 Concluding remarks 

We have used SPH and the software package mmas to study triple- 
star collisions. Although such collisions span a tremendous amount 
of parameter space, our modest number of SPH calculations do pro- 
vide some valuable insights. For the (parabolic) encounters that we 
consider, we find that the order in which stars collide (sec §3.1.1), 
the angle of approach of the third star (§3.1.2), and the periastron 
separation of the collisions (§3.1.3) have only a slight effect on the 
chemical composition distribution within the final collision prod- 
uct. The order and orbital parameters of the collisions can, however, 
significantly affect the size and structure of the product. 

The results of §3.2.1 help establish that the simple fluid sort- 
ing algorithm of mmas reproduce the important features of our SPH 
models, even when one of the parent stars is itself a collision prod- 
uct. The MMAS package can therefore be considered an adequate, 
if not an accurate, substitute for a hydrodynamics code in many 
situations. This realization will help simplify the process of gen- 
erating collision product models in cluster simulations, because a 
full hydrodynamics calculation will not necessarily need to be run 
for each collision. Indeed, we hope the mmas package will be used 
to help account for stellar collisions in dynamics simulations of 
globular clusters. Toward this end, MMAS is already being incorpo- 
rated into two software packages, TRIPTYCH and tripletych, that 
respectively treat encounters between two stars and among three 
stars (see Sills et al. 2003).^ These packages are controlled through 
a web interface and treat the orbital trajectories, possible merger(s), 
and evolution of the merger product and therefore incorporate three 
main branches of stellar astrophysics: dynamics, hydrodynamics, 
and evolution. 

^ See http://faculty.vassar.edu/lorribardi/triptych/ 
and http: //f acuity .vassar. edu/lombardi /tripletych/. 



The product size estimates of §3.2.2 are admittedly crude. For 
example, partial ionization and radiation pressure are neglected. Al- 
though the exact size of a collision product is difficult to determine, 
our calculations indicate that the first and final collision products 
are always significantly larger than their thermally relaxed coun- 
terparts would be. Indeed, according to our MMAS calculations in 
§3.2.2, the final collision product can have a radius up to ^ lO^R©, 
easily exceeding the size of a typical red giant. Furthermore, these 
calculations have assumed that some mechanism has braked the of- 
ten rapid rotation of the products, so any rotation that does remain 
will only further enhance the size of the products. The extended 
sizes of the products will increase the multi-star collision rate over 
that calculated in previous treatments of binary-single and binary- 
binary encounters. 

All of the scenarios we consider with SPH in this paper in- 
volve one 0.8 and two O.6M0 stars. Without shock heating, the 
low-A fluid of the O.SM© star would sink to the core of the final 
collision product, while its high-^ portions of the O.SMq would 
settle in the outer layers. The intermediate layers of the product 
would consist of fluid with the same A range from all the parents. 
Simply sorting the fluid in this way, without running a hydrody- 
namics calculation, can therefore provide a zeroth-order model of 
the collision product that captures some of its important qualitative 
features (see Fig. 3). 

However, non-uniform shock heating during the collisions 
somewhat alters the relative values of the entropic variable A in the 
fluid, resulting in a slightly different sorting pattern. Because the 
amount and distribution of shock heating are dependent on the de- 
tails of a collision, the sorting of the fluid varies with, for example, 
the order in which the stars collide. Shock heating can have larger 
consequences on the chemical abundance profiles of elements, such 
as C^^, that exist in substantial amounts only in a small shell in 
the initial parent stars. However, the chemical abundance profiles 
of most elements, particularly helium, are always qualitatively the 
same, regardless of how the three stars are merged. Because the 
abundance and distribution of helium (and hence hydrogen) is one 
of the most important factors in determining the collision product's 
subsequent course of stellar evolution, we believe that the order 
and geometry of the collisions will not significantly affect the stel- 
lar evolution of the product. Indeed, Sills et al. (2002) have recently 
presented a set of stellar evolution calculations for a collision prod- 
uct for which the starting yrec models were generated from SPH 
calculations of different resolutions. The variations in the initial he- 
lium profiles of their models are roughly comparable to those in our 
helium profiles resulting from colliding three parent stars in differ- 
ent ways. Although Sills et al. (2002) do find detailed differences in 
the evolution, especially during the "pre-main-sequence" contrac- 
tion, the evolutionary tracks and time- scales are quite similar. We 
therefore feel that, for low-velocity collisions, the hydrodynamical 
details of how three stars are merged will not significantly affect 
the stellar evolution of the collision product — the major caveat here 
being that the geometry of the collisions can of course affect the ro- 
tation of the product, which in turn can greatly affect its evolution 
(Sills et al. 2001). 

Surface abundances of lithium and beryllium are particularly 
interesting to monitor, as these elements can be used as observa- 
tional indicators of mixing and perhaps coUisional history. As in the 
single-single star collisions presented by Lombardi et al. (2002), we 
find that the triple-star mergers presented here yield collision prod- 
ucts that are severely depleted of lithium and somewhat of beryl- 
lium at the surface. Even in the relatively gentle (parabolic) cases 
that we have considered, the collisions are energetic enough to ex- 
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pel most of the lithium and beryllium from the outer layers of the 
parents. 

4.2 Future work 

There are many scenarios to explore when dealing with collisions in 
environments as chaotic as dense stellar systems. Different orbital 
geometries besides the parabolic trajectories treated in this paper 
still need to be considered in more detail. Large stellar velocities 
in galactic nuclei lead to hyperbolic collisions. In globular clusters, 
perturbations to a binary can lead to an elliptical collision, while an 
encounter with a very hard binary can lead to significantly hyper- 
bolic collisions. 

Future studies may want to include a more detailed look at the 
hydrodynamics during grazing encounters, which could be done 
efficiently with the help of GRAPE (short for GRAvity PipE) spe- 
cial purpose hardware for calculating the self-gravity of the system. 
Furthermore, encounters involving more than three stars, such as in 
binary-binary interactions, may also warrant further examination: 
for example, the final collision product generated in a triple-star 
merger is typically so extended that it could immediately start suf- 
fering Roche lobe overflow if left in orbit around a fourth star 

Collisions among a larger variety of stellar types and masses, 
reflective of the diverse populations of clusters, will also need 
to be explored. We have been concentrating on low-mass main- 
sequence stars, but collisions between high-mass main-sequence 
stars in young compact star clusters, or giants located in the dense 
cores of globular clusters, for example, are frequent. A logical first 
step would be to examine high-mass main- sequence stars in a run- 
away merger scenario. It would therefore be very useful to develop 
a generalization of the fluid sorting method that includes radiation 
pressure in the equation of state. 

Due to shock heating during the collision, the product is much 
larger than a thermally equilibrated main-sequence star of the same 
mass. How much of an effect this increased radius has on the effec- 
tive cross section for merger is subject to many variables, includ- 
ing the structure of the product's outer layers and the velocity of 
approaching stars. In environments such as active galactic nuclei, 
where relative velocities tend to be high, the low-density outer lay- 
ers of a newly formed collision product could Ukely get stripped 
by passing stars. However, in globular clusters, where stellar veloc- 
ities tend to be small, collisions with even low-density envelopes 
may lead to significantly increased rates of merger. It would be use- 
ful to develop a robust collision module that could quickly predict 
whether any given collision trajectory will lead to a merger, and, if 
not, describe how the stars are affected by the interaction. 

One simple approximation often implemented in cluster sim- 
ulations is that a collision product instantaneously achieves its ther- 
mally relaxed radius, a good approximation when the time between 
collisions is much longer than the thermal time-scale. Arguing in- 
stead that the global thermal time-scale of the first product can be 
much larger than the time between collisions in interactions in- 
volving binaries, we make a different approximation in this pa- 
per, namely that the first product's radius (and more generally its 
structure) does not substantially evolve between collisions. Future 
scattering experiments could model thermally relaxing stars and 
study more carefully the timescale between collisions mediated by 
binaries. The thermal time-scale in the outer layers of a collision 
product can be orders of magnitude less than its global thermal 
time-scale (see table 1 of Sills et al. 1997), so that it may actu- 
ally be necessary to follow the thermal contraction and stellar or- 
bits simultaneously. Indeed, in the extremely low density layers of 



a collision product, it is even possible for the thermal time-scale 
to be comparable to the (hydro)dynamical time-scale, so that the 
product could undergo significant thermal contraction even before 
it reaches hydrodynamical equilibriiun. It would be helpful if fu- 
ture stellar evolution calculations of collision products included a 
detailed description of the products' size and structure throughout 
the thermal relaxation stage. How quickly the outer layers of the 
thermally expanded product change with time will substantially af- 
fect its likeUhood of subsequent collisions. Initial conditions for 
such stellar evolution calculations could be provided by the pub- 
licly available mmas package. 

The primary hurdle for incorporating collisions into realistic 
stellar dynamics simulations is currently the stellar evolution of the 
collision products. Such stars are highly non-canonical, typically 
with very peculiar structural and composition profiles, and present 
a challenging set of initial conditions for stellar evolution codes. 
To make matters even more intricate, rotation, which is typically 
rapid after merging, will affect the structural properties and chemi- 
cal compositions of the stars as they evolve (e.g.. Sills et al. 2001). 
This rapid rotation also has the effect of ejecting mass as the prod- 
uct thermally contracts. Studying this emitted mass will be worth- 
while, as it may likely carry away angular momentum and at least 
partially brake rotating collision products. 
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